On cherche à étudier l’effet de trois facteurs sur le transcriptome des racines d’Arabidopsis thaliana. Le CO2, au cours des études préliminaires, s’est montré peu influent en conditions contrôles de fer et de nitrates, et accentué en cas de stress nutritionnel. Nous reprennons ces résultats avec des fonctions génériques et propres pour en faire le résumé et de jolis graphes.
Import des données : matrice d’expression
On a, pour chaque gène et chaque condition, son niveau d’expression en sortie de quantification. On labelle les conditions avec le code suivant : lettre majuscule pour le niveau fort, minuscule pour le niveau faible. Le réplicat est donné après l’underscore.
# quantification file
data <- read.csv("quantifFiles/QuantifGenes.csv", h = T, sep = ",")
rownames(data) <- data$Gene
genes = which(!(grepl("__", rownames(data))))
not_quant = data[which((grepl("__", rownames(data)))), ]
data = data[genes, grepl("R", colnames(data))]
getLabel("R6")[1] "cnF_3"[1] At_AmbientCO2_LowNitrate_Fe1
48 Levels: At_AmbientCO2_HighNitrate_Fe1 ... Sl_ElevatedCO2_LowNitrate_FeStarvation3[1] "At_AmbientCO2_LowNitrate_Fe"keep <- rowSums(data) >= 10
data <- data[keep, ]
group <- sapply(colnames(data), getLabel, with.rep = F)
colnames(data) <- sapply(colnames(data), getLabel)
head(data) cNF_3 cNF_2 cNF_1 cnF_2 cnF_1 cnF_3 CNF_1 CnF_2 CnF_1 CnF_3 cNf_1
AT1G01010 1526 1006 1116 1275 967 1018 854 1132 1294 1364 2325
AT1G01020 416 285 289 349 364 370 226 513 502 561 461
AT1G01030 31 15 19 29 36 28 12 47 34 47 18
cnf_2 cnf_1 cNf_2 cNf_3 cnf_3 Cnf_3 CNf_1 Cnf_1 Cnf_2 CNf_3 CNf_2
AT1G01010 2113 2193 2564 2364 2074 1987 2027 1754 1697 1537 1898
AT1G01020 407 432 614 380 502 484 467 426 415 413 462
AT1G01030 40 32 44 37 27 42 39 36 40 37 37
CNF_3 CNF_2
AT1G01010 816 912
AT1G01020 223 312
AT1G01030 15 19
[ reached 'max' / getOption("max.print") -- omitted 3 rows ][1] 23342 24Effet CO2 en conditions contrôle
On définit les conditions contrôle comme suit : fort nitrate et fort fer.
method = "edger"
g = list()
# reference condition as first element of labels
labels <- c("cNF", "CNF")
genes1 <- dualDE(data, labels, pval = 0.01, method = method) (Intercept) groupCNF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cNF_3" "cNF_2" "cNF_1" "CNF_1" "CNF_3" "CNF_2"
cNF_3 cNF_2 cNF_1 CNF_1 CNF_3 CNF_2
1.0090325 1.0047763 0.9805842 1.0116521 0.9699686 1.0239863 
[1] "68 genes DE"
[1] 17Low nitrate
(Intercept) groupCnF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cnF_2" "cnF_1" "cnF_3" "CnF_2" "CnF_1" "CnF_3"
cnF_2 cnF_1 cnF_3 CnF_2 CnF_1 CnF_3
0.9862253 0.9873030 0.9924054 1.0167225 1.0047221 1.0126218 
[1] "1312 genes DE"
Low iron
(Intercept) groupCNf
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cNf_1" "cNf_2" "cNf_3" "CNf_1" "CNf_3" "CNf_2"
cNf_1 cNf_2 cNf_3 CNf_1 CNf_3 CNf_2
0.9799002 1.0044883 0.9723876 0.9981139 1.0241404 1.0209697 
[1] "2219 genes DE"
Low iron and nitrate
(Intercept) groupCnf
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cnf_2" "cnf_1" "cnf_3" "Cnf_3" "Cnf_1" "Cnf_2"
cnf_2 cnf_1 cnf_3 Cnf_3 Cnf_1 Cnf_2
0.9975633 1.0046643 1.0015093 0.9727479 1.0064673 1.0170479 
[1] "1656 genes DE"
Venn diagram
On visualise les gènes différentiellement exprimés en commun entre les différents niveaux des autres facteurs.
library(ggVennDiagram)
library(VennDiagram)
gene_list <- list()
for (comp in names(g)) {
print(g[[comp]])
gene_list[[comp]] <- g[[comp]]$gene_id
} gene_id a.value m.value p.value q.value rank estimatedDEG
1 AT5G64120 8.797218 -1.1121447 1.895001e-15 4.423312e-11 1 1
2 AT1G19250 9.948535 -0.8706103 1.978749e-14 1.945820e-10 2 1
3 AT4G25100 9.324940 -0.9038503 2.500840e-14 1.945820e-10 3 1
4 AT5G59520 8.819711 -1.2066060 1.692883e-12 9.014405e-09 4 1
5 AT2G01520 14.038784 -0.9716064 1.930941e-12 9.014405e-09 5 1
6 AT1G67270 5.104540 -2.2140173 3.322466e-12 1.292550e-08 6 1
7 AT5G23980 8.719121 -1.1156310 2.014654e-10 6.718007e-07 7 1
8 AT4G37060 10.229806 -0.8191550 6.994823e-10 2.040914e-06 8 1
9 AT5G23990 4.823373 -2.0335514 9.538085e-10 2.473755e-06 9 1
upreg
1 0
2 0
3 0
4 0
5 0
6 0
7 0
8 0
9 0
[ reached 'max' / getOption("max.print") -- omitted 59 rows ]
gene_id a.value m.value p.value q.value rank estimatedDEG
1 AT5G61600 10.978022 1.715345 5.773443e-56 1.347637e-51 1 1
2 AT4G27450 10.987889 1.361106 1.242962e-51 1.450661e-47 2 1
3 AT1G27730 9.421760 2.187544 2.461130e-51 1.914923e-47 3 1
4 AT5G51190 8.417337 2.782680 4.320191e-49 2.521047e-45 4 1
5 AT2G27830 9.879502 1.316552 1.397616e-39 6.524629e-36 5 1
6 AT1G19530 12.542251 1.538979 1.886217e-38 7.338012e-35 6 1
7 AT4G27280 10.488397 1.691216 4.611327e-38 1.537680e-34 7 1
8 AT5G13930 10.407278 -1.747350 6.722757e-38 1.961533e-34 8 1
9 AT3G44260 8.737846 2.483672 1.819583e-36 4.719190e-33 9 1
upreg
1 1
2 1
3 1
4 1
5 1
6 1
7 1
8 0
9 1
[ reached 'max' / getOption("max.print") -- omitted 1303 rows ]
gene_id a.value m.value p.value q.value rank estimatedDEG
1 AT1G78340 10.002535 -3.261719 5.399049e-136 1.260246e-131 1 1
2 AT4G13180 10.301728 -2.811270 3.011796e-121 3.515067e-117 2 1
3 AT4G15550 8.217350 -3.781600 5.522983e-117 4.297249e-113 3 1
4 AT4G34131 8.946763 -2.699004 1.584031e-95 9.243612e-92 4 1
5 AT1G76680 9.778808 -3.051667 2.789158e-90 1.302091e-86 5 1
6 AT3G28345 9.654586 -2.538451 3.170729e-89 1.233519e-85 6 1
7 AT3G45060 9.161753 2.510301 5.415418e-82 1.805810e-78 7 1
8 AT1G76690 9.762616 -2.235867 5.461401e-76 1.593500e-72 8 1
9 AT3G28740 4.011621 -7.093650 9.075419e-76 2.353760e-72 9 1
upreg
1 0
2 0
3 0
4 0
5 0
6 0
7 1
8 0
9 0
[ reached 'max' / getOption("max.print") -- omitted 2210 rows ]
gene_id a.value m.value p.value q.value rank estimatedDEG
1 AT3G49620 8.971142 3.325629 4.576332e-77 1.068207e-72 1 1
2 AT2G26400 12.032573 -2.107176 1.946526e-64 2.271791e-60 2 1
3 AT1G21100 9.986600 1.996729 3.562033e-43 2.771499e-39 3 1
4 AT1G69530 9.527594 -1.507353 2.066807e-40 1.206085e-36 4 1
5 AT5G51860 6.042441 3.573607 4.906245e-40 2.290431e-36 5 1
6 AT5G51870 8.032896 1.826646 6.948595e-39 2.703235e-35 6 1
7 AT1G08630 8.692334 2.049414 1.918515e-37 6.397426e-34 7 1
8 AT1G56600 8.952644 -1.850385 1.108793e-36 3.235181e-33 8 1
9 AT1G43800 11.887403 -2.934754 2.086789e-36 5.412203e-33 9 1
upreg
1 1
2 0
3 1
4 0
5 1
6 1
7 1
8 0
9 0
[ reached 'max' / getOption("max.print") -- omitted 1647 rows ]
partitions <- get.venn.partitions(gene_list, keep.elements = T)
partitions$shared <- rowSums(partitions[1:4])
partitions <- partitions[order(-partitions$shared), ]
# common_genes <- unlist(partitions[1, '..values..']) results <- getBM( filters =
# 'ensembl_gene_id', attributes = c('ensembl_gene_id', 'description',
# 'external_gene_name', 'entrezgene_id'), values = common_genes, mart = mart)
# results <- results[!rownames(results) %in%
# which(duplicated(results$ensembl_gene_id)), ] kable(results)
# sharedBy4 <- unique(unlist(subset(partitions,
# partitions$shared==3)$..values..)) a <- OntologyProfileAt(sharedBy4)Summary graph
# getExpression('AT2G35980', conds = c('cNF', 'CNF'))
d <- data.frame(matrix(ncol = 4, nrow = 2))
colnames(d) <- names(g)
row.names(d) <- c("up", "down")
for (comp in names(g)) {
d["up", comp] <- sum(g[[comp]]$upreg == 1)
d["down", comp] <- sum(g[[comp]]$upreg == 0)
}
res <- melt(d)
res$reg = rep(c("up", "down"), 4)
genesco2At <- res
save(genesco2At, file = "genesco2At.RData")
ggplot(res, aes(fill = reg, y = value, x = variable)) + geom_bar(position = "stack",
stat = "identity", alpha = 0.5, color = "black") + ggtitle("Carbon dioxide effet on gene regulation") +
xlab("") + ylab("Number of differentially expressed genes") + coord_flip()
Investigations
Double carence en CO2 ambient et élevé
On s’interroge sur l’effet de la double carence fer et nitrate en CO2 ambiant et élevé. On ne l’avait pas fait car cela fait varer deux facteurs à la fois.
(Intercept) groupcNF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cnf_2" "cnf_1" "cnf_3" "cNF_3" "cNF_2" "cNF_1"
cnf_2 cnf_1 cnf_3 cNF_3 cNF_2 cNF_1
1.0242524 1.0393358 1.0253221 0.9831632 0.9914864 0.9364401 
[1] "9069 genes DE"
(Intercept) groupCNF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "Cnf_3" "Cnf_1" "Cnf_2" "CNF_1" "CNF_3" "CNF_2"
Cnf_3 Cnf_1 Cnf_2 CNF_1 CNF_3 CNF_2
0.9823882 1.0318759 1.0281336 1.0008911 0.9447685 1.0119428 
[1] "8309 genes DE"
doubleInt <- list(`ambiant CO2` = genes5$gene_id, `elevanted CO2` = genes6$gene_id)
ggVennDiagram(doubleInt)
ensembl_gene_id
1 AT1G01720
2 AT1G01750
3 AT1G02575
4 AT1G02640
5 AT1G02850
6 AT1G03790
7 AT1G03870
8 AT1G04250
9 AT1G05650
10 AT1G05660
11 AT1G06120
12 AT1G07680
13 AT1G07690
14 AT1G07750
15 AT1G08080
16 AT1G08100
17 AT1G08290
18 AT1G08500
description
1 NAC domain-containing protein 2 [Source:UniProtKB/Swiss-Prot;Acc:Q39013]
2 actin depolymerizing factor 11 [Source:TAIR;Acc:AT1G01750]
3 unknown protein; BEST Arabidopsis thaliana protein match is: unknown protein (TAIR:AT1G02570.1). [Source:TAIR;Acc:AT1G02575]
4 Probable beta-D-xylosidase 2 [Source:UniProtKB/Swiss-Prot;Acc:Q94KD8]
5 Beta-glucosidase 11 [Source:UniProtKB/Swiss-Prot;Acc:B3H5Q1]
6 Zinc finger CCCH domain-containing protein 2 [Source:UniProtKB/Swiss-Prot;Acc:Q9ZWA1]
7 At1g03870 [Source:UniProtKB/TrEMBL;Acc:B3LF88]
8 Auxin-responsive protein IAA17 [Source:UniProtKB/Swiss-Prot;Acc:P93830]
9 F3F20.10 protein [Source:UniProtKB/TrEMBL;Acc:Q9SYK6]
10 F3F20.11 protein [Source:UniProtKB/TrEMBL;Acc:Q9SYK7]
11 Delta-9 desaturase-like 3 protein [Source:UniProtKB/Swiss-Prot;Acc:Q9FPD5]
12 F24B9.23 [Source:UniProtKB/TrEMBL;Acc:Q9LQP3]
13 FUNCTIONS IN: molecular_function unknown; INVOLVED IN: biological_process unknown; LOCATED IN: endomembrane system; EXPRESSED IN: root; CONTAINS InterPro DOMAIN/s: Orthoreovirus membrane fusion p10 (InterPro:IPR009854); BEST Arabidopsis thaliana pro /.../atch is: unknown protein (TAIR:AT1G54950.1); Ha. [Source:TAIR;Acc:AT1G07690]
14 At1g07750/F24B9_13 [Source:UniProtKB/TrEMBL;Acc:Q9LQQ3]
15 Alpha carbonic anhydrase 7 [Source:UniProtKB/Swiss-Prot;Acc:Q8L817]
16 NRT2 [Source:UniProtKB/TrEMBL;Acc:A0A178W7F7]
17 Zinc finger protein WIP3 [Source:UniProtKB/Swiss-Prot;Acc:Q9SGD1]
18 Early nodulin-like protein 18 [Source:UniProtKB/TrEMBL;Acc:O82083]
external_gene_name entrezgene_id
1 NAC002 839265
2 ADF11 839281
3 10723100
4 BXL2 837940
5 BGLU11 839435
6 SOM 839408
7 FLA9 839384
8 IAA17 839568
9 837072
10 837073
11 837121
12 837281
13 837282
14 837290
15 ACA7 837326
16 NRT2.2 837328
17 WIP3 837349
18 ENODL18 837371
[ reached 'max' / getOption("max.print") -- omitted 6632 rows ]On trouve énormément de gènes, ce qui est cohérent car le nitrate et le fer pris séparément avaient déjà beaucoup d’effet.
DEG cNF et CNF avec un autre condition
On cherche à savoir si il n’y pas pas un soucis avec l’échantillon CNF, qui semble avoi très très peu de différences avec cNF. Aurait-on envoyé les mêmes tubes?
On ne voit que très peu de DEGs en comparant directement cNF et CNF, on cherche donc à savoir si on trouve le même nombre de DEGs quand on compare cNF - x et CNF -x, ce qui validerait l’hypothèse d’un transcriptôme quasi identique. On choisit ici par exemple \(x=\) cnF.
(Intercept) groupcnF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cNF_3" "cNF_2" "cNF_1" "cnF_2" "cnF_1" "cnF_3"
cNF_3 cNF_2 cNF_1 cnF_2 cnF_1 cnF_3
0.9859482 0.9865531 0.9545113 1.0150637 1.0299591 1.0279646 
[1] "2803 genes DE"
(Intercept) groupcnF
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "CNF_1" "CNF_3" "CNF_2" "cnF_2" "cnF_1" "cnF_3"
CNF_1 CNF_3 CNF_2 cnF_2 cnF_1 cnF_3
0.9885185 0.9468228 1.0077432 1.0150001 1.0207525 1.0211628 
[1] "1928 genes DE"
double <- list(`ambiant CO2` = genes7$gene_id, `elevanted CO2` = genes8$gene_id)
ggVennDiagram(double)
Il semble qu’il y ait quand même pas mal de différences entre ces transcriptomes, car quand on les compare au même troisième transcriptome on a une différence de 1000 gènes (1900 contre 2800). Le diagramme de Venn montre que ces gènes sont en partie différents. Les transcriptômes ont bien l’air différents. Si on prend un autre x :
(Intercept) groupcNf
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "cNF_3" "cNF_2" "cNF_1" "cNf_1" "cNf_2" "cNf_3"
cNF_3 cNF_2 cNF_1 cNf_1 cNf_2 cNf_3
0.9844303 0.9894266 0.9466501 1.0201203 1.0475235 1.0118491 
[1] "7341 genes DE"
(Intercept) groupcNf
1 1 0
2 1 0
3 1 0
4 1 1
5 1 1
6 1 1
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"
[1] "CNF_1" "CNF_3" "CNF_2" "cNf_1" "cNf_2" "cNf_3"
CNF_1 CNF_3 CNF_2 cNf_1 cNf_2 cNf_3
0.9891655 0.9317337 1.0004327 1.0190222 1.0502432 1.0094027 
[1] "7140 genes DE"
double <- list(`ambiant CO2` = genes7$gene_id, `elevanted CO2` = genes8$gene_id)
ggVennDiagram(double)
ensembl_gene_id
1 AT1G01190
2 AT1G01720
3 AT1G02575
4 AT1G02850
5 AT1G03870
6 AT1G05660
7 AT1G05680
8 AT1G06120
9 AT1G07680
10 AT1G07690
11 AT1G09750
12 AT1G09932
13 AT1G10385
14 AT1G11545
15 AT1G12030
16 AT1G12080
17 AT1G12090
18 AT1G12200
description
1 CYP78A8 [Source:UniProtKB/TrEMBL;Acc:A0A178WIC4]
2 NAC domain-containing protein 2 [Source:UniProtKB/Swiss-Prot;Acc:Q39013]
3 unknown protein; BEST Arabidopsis thaliana protein match is: unknown protein (TAIR:AT1G02570.1). [Source:TAIR;Acc:AT1G02575]
4 Beta-glucosidase 11 [Source:UniProtKB/Swiss-Prot;Acc:B3H5Q1]
5 At1g03870 [Source:UniProtKB/TrEMBL;Acc:B3LF88]
6 F3F20.11 protein [Source:UniProtKB/TrEMBL;Acc:Q9SYK7]
7 Glycosyltransferase (Fragment) [Source:UniProtKB/TrEMBL;Acc:A0A0K1SBE8]
8 Delta-9 desaturase-like 3 protein [Source:UniProtKB/Swiss-Prot;Acc:Q9FPD5]
9 F24B9.23 [Source:UniProtKB/TrEMBL;Acc:Q9LQP3]
10 FUNCTIONS IN: molecular_function unknown; INVOLVED IN: biological_process unknown; LOCATED IN: endomembrane system; EXPRESSED IN: root; CONTAINS InterPro DOMAIN/s: Orthoreovirus membrane fusion p10 (InterPro:IPR009854); BEST Arabidopsis thaliana pro /.../atch is: unknown protein (TAIR:AT1G54950.1); Ha. [Source:TAIR;Acc:AT1G07690]
11 Aspartyl protease AED3 [Source:UniProtKB/Swiss-Prot;Acc:O04496]
12 Phosphoglycerate mutase family protein [Source:UniProtKB/TrEMBL;Acc:Q8GWG7]
13 Exocyst complex component EXO84A [Source:UniProtKB/Swiss-Prot;Acc:F4I4B6]
14 Xyloglucan endotransglucosylase/hydrolase [Source:UniProtKB/TrEMBL;Acc:A0A178W0W2]
15 At1g12030 [Source:UniProtKB/TrEMBL;Acc:O65376]
16 At1g12080 [Source:UniProtKB/TrEMBL;Acc:O65370]
17 ELP [Source:UniProtKB/TrEMBL;Acc:A0A178W1A4]
18 Flavin-containing monooxygenase [Source:UniProtKB/TrEMBL;Acc:A0A178W1K3]
external_gene_name entrezgene_id
1 CYP78A8 839233
2 NAC002 839265
3 10723100
4 BGLU11 839435
5 FLA9 839384
6 837073
7 UGT74E2 837075
8 837121
9 837281
10 837282
11 AED3 837504
12 837526
13 EXO84A 837578
14 XTH8 837698
15 837755
16 837760
17 ELP 837761
18 837773
[ reached 'max' / getOption("max.print") -- omitted 6047 rows ]On va faire les corrélations d’expression entre ces deux conditions.
scatter <- function(df, x, y) {
sp <- ggplot(df, aes(log(df[, x]), log(df[, y]))) + geom_bin2d(bins = 70) + scale_fill_gradient2() +
theme(legend.position = "none") + labs(x = x, y = y)
return(sp)
}
comps <- labels <- c("cNF", "CNF")
diffs <- c()
# correlations between ambient and elevated
for (ambient in colnames(data)[grepl("cNF", colnames(data))]) {
for (elevated in colnames(data)[grepl("CNF", colnames(data))]) {
diffs = c(diffs, cor(data[ambient], data[elevated]))
}
}
diffs_other <- c()
for (ambient in colnames(data)[grepl("cNF", colnames(data))]) {
for (elevated in colnames(data)[!grepl("NF", colnames(data))]) {
diffs_other = c(diffs_other, cor(data[ambient], data[elevated]))
}
}
# correlations within ambien and elevated
sames <- c()
samples <- colnames(data)[grepl("cNF", colnames(data))]
sames <- c(sames, cor(data[, samples[1]], data[, samples[2]]), cor(data[, samples[1]],
data[, samples[3]]), cor(data[, samples[2]], data[, samples[3]]))
samples <- colnames(data)[grepl("CNF", colnames(data))]
sames <- c(sames, cor(data[, samples[1]], data[, samples[2]]), cor(data[, samples[1]],
data[, samples[3]]), cor(data[, samples[2]], data[, samples[3]]))
sames[1] 0.9938199 0.9906937 0.9810858 0.9860430 0.9921428 0.9860068df <- data.frame(`Same condition` = c(sames, rep(NA, 3)), `Ambient vs Elevated` = diffs)
res <- na.omit(melt(df))
# comparisons
ggplot(data = res, aes(x = variable, y = value, fill = variable)) + geom_boxplot(binaxis = "y",
alpha = 0.4, stackdir = "center") + theme(axis.text.x = element_text(angle = 320,
hjust = 0, colour = "grey50"), plot.title = element_text(size = 14, face = "bold")) +
ggtitle("Correlation for expression accross conditions") + geom_dotplot(binaxis = "y",
stackdir = "center") + stat_compare_means(method = "wilcox.test", hide.ns = FALSE,
label = "p.signif", symnum.args = list(cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05,
1), symbols = c("****", "***", "**", "*", "ns")))
res <- rbind.data.frame(res, data.frame(variable = rep("cNF vs all others except CNF",
length(diffs_other)), value = diffs_other))
ggplot(data = res, aes(x = variable, y = value, fill = variable)) + geom_boxplot(binaxis = "y",
alpha = 0.4, stackdir = "center") + theme(axis.text.x = element_text(angle = 320,
hjust = 0, colour = "grey50"), plot.title = element_text(size = 14, face = "bold")) +
ggtitle("Correlation for expression accross conditions") + geom_dotplot(binaxis = "y",
stackdir = "center") + stat_compare_means(method = "wilcox.test", hide.ns = FALSE,
label = "p.signif", symnum.args = list(cutpoints = c(0, 1e-04, 0.001, 0.01, 0.05,
1), symbols = c("****", "***", "**", "*", "ns")))
cNF et CNF semblent bien plus proches entre eux que cNF ne l’est avec tous les autres.